par(pch=20)
n=5000
x_true_sd=10; x_true_mean=50
y_true_sd=10; y_true_mean=60
X=rnorm(n,mean=x_true_mean,sd=x_true_sd)
Y=rnorm(n,mean=y_true_mean,sd=y_true_sd)

#---
# *** populations
plot(c(X,Y),col=c(rep(3,n),rep(4,n)))
abline(h=mean(X),col=3)
abline(h=x_true_mean,col=3,lty=2)
abline(h=x_true_mean+c(-1,1)*x_true_sd,col=3,lty=3)
abline(h=mean(Y),col=4)
abline(h=y_true_mean,col=4,lty=2)
abline(h=y_true_mean+c(-1,1)*y_true_sd/2,col=4,lty=3)
#----
# *** single run showing estimated parameters vs. sample size
# *** 
samples=c()
for (sample_size in seq(1,n,20)) {
	x=sample(X,size=sample_size,replace=F)
	sample_mean=mean(x)
	sample_sd=sd(x)
	samples=rbind(samples,
		c(size=sample_size,mean=sample_mean,sd=sample_sd))
}
samples=data.frame(samples)
par(mfrow=c(2,1))
plot(((mean-x_true_mean)/x_true_mean)~size,data=samples,
	xlab='sample size',ylab='% deviation from true mean',
	ylim=c(-1,1)/2,main="sample mean vs. population mean")
abline(h=0,col=3)
plot(((sd-x_true_sd)/x_true_sd)~size,data=samples,
	xlab='sample size',ylab='% deviation from true sd',
	ylim=c(-1,1)/2,main="sample sd vs. population sd")
abline(h=0,col=3)
#----
# create artificial variates
beta_x0=7.621
beta_x1=1.973
beta_y1=9.648
beta_x2=0.123

Z=beta_x0+beta_x1*X+beta_x2*X+beta_y1*Y
d=data.frame(cbind(X=X,Y=Y,Z=Z))

#----
#x=c(8.6, 9.4, 7.9, 6.8, 8.3, 7.3, 9.2, 9.6, 8.7, 11.4, 10.3, 5.4, 8.1, 5.5, 6.9)
sim=function (n,mean,sd) {
	#x=rnorm(n=n,mean=mean,sd=sd)
	x=runif(n=n,min=mean-3*sd,max=mean+3*sd)
	r=c(); for (i in 1:length(x)) { 
		r=rbind(r,
			c(
				n=i,
				mu=mean(x[1:i]),
				sigma=sd(x[1:i]),
				se=sd(x[1:i])/sqrt(i)
			)
		) 
	}
	r=data.frame(r)
	r=cbind(r,mu_err=r$mu-mean(x),sigma_err=r$sigma-sd(x))
	r=cbind(r,se=r$sigma/sqrt(r$n))
	r=cbind(r,err=r$se * qt(0.975,df=r$n-1))
	plot.data=cbind(x, mu=r$mu, 
		low.mu=r$mu-r$err, hi.mu=r$mu+r$err,
		low.mu.s=r$mu-r$sigma,hi.mu.s=r$mu+r$sigma)
	col=c(2,1,1,1,3,3)
	typ=c('p',rep('l',5))
	lty=rep(1,5)
	matplot(plot.data,col=col,typ=typ,lty=lty,pch=20)
	abline(h=c(mean(x)-sd(x),mean(x),mean(x)+sd(x)),lty=c(1,2,1),col=2)
	legend(length(x)/2,min(x)-2*sd(x),colnames(plot.data),col=col,lty=lty,cex=0.7)
}

# *** Exchange rates
d=read.table(file='http://www.stat.duke.edu/data-sets/mw/ts_data/spot_exrates.dat')
colnames(d)=c('AUD','BEF','CAD','FRF','DEM','JPY','NLG','NZD','ESP','SEK','CHF','GBP')
d=ts(d)
pairs(d,pch='.',col=4,main='Original series')
dd1=diff(d)
pairs(dd1,pch='.',col=4,main='Differenced series (t-1)')
plot(d[,1:5])
plot(cbind(GBP=d[,'GBP'],GBPd1=dd1[,'GBP']),main='GBP original and t-1 series')

# -- cut data set to GBP
f <- function(symbol) { 
	dx=data.frame(v=d[,symbol])
	dx=cbind(t=row(dx),dx) # add row number
	par(mfrow=c(3,3))
	for (i in c(50,100,200,500,1000,1500,2000,2200,2500)) {
		dx.sample=dx[1:i,] 
		fit=lm(v~t,data=dx.sample)
		dx.pred=predict(fit,data.frame(t=(i+1):nrow(dx)))
		data=cbind(dx,fitted=c(fit$fitted.values,dx.pred))
		matplot(data[,2:3],typ='l',ylab='',lty=1,col=3:4)
		abline(v=i,col=4)
		#readline()
	}
}

# *** Artificial ARIMA
show.arima <- function(data,order,p) {
	train.samples=floor(length(data) * p)
	fit=arima(data[1:train.samples],order=order)
	pred=predict(fit,n.ahead=length(data)-train.samples)$pred
	error=data[(train.samples+1):length(data)]-pred
	pad=rep(NA,train.samples)
	matplot(
		cbind(
			data,
			c(pad,pred),
			c(pad,error)
		),
		typ='l',pch=20,cex=0.3,cex.main=1,ylab=NA,
		main=sprintf("(p=%d,d=%d,q=%d)\nRMSE=%3.2f , sigma^2=%3.2f\nloglik=%3.2f , aic=%3.2f",
			order[1],order[2],order[3],sqrt(sum(fit$resid^2/25)),fit$sigma2,fit$loglik,fit$aic
		)
	)
	abline(v=train.samples,col=2)
	return(fit)
}

graphics.off()
x=1:50+rep(c(1,2,3,2,-1,2,1,4,2,1),5)
x=jitter(x,amount=1)
windows()
par(mfrow=c(4,1))
plot(x,typ='b',pch=20,cex=0.5)
plot(diff(x),typ='b',pch=20,cex=0.5)
acf(x)
pacf(x)
windows()
par(mfrow=c(3,3),mai=c(0.1,0.1,0.1,0.1)*4)
p=0.5
fits=list()
fits=c(fits,list(show.arima(x,order=c(p=0,d=0,q=0),p=p)))
fits=c(fits,list(show.arima(x,order=c(p=0,d=1,q=1),p=p)))
fits=c(fits,list(show.arima(x,order=c(p=0,d=2,q=2),p=p)))
fits=c(fits,list(show.arima(x,order=c(p=1,d=2,q=2),p=p)))
fits=c(fits,list(show.arima(x,order=c(p=2,d=2,q=2),p=p)))
fits=c(fits,list(show.arima(x,order=c(p=4,d=2,q=2),p=p)))
fits=c(fits,list(show.arima(x,order=c(p=5,d=2,q=3),p=p)))
fits=c(fits,list(show.arima(x,order=c(p=5,d=3,q=3),p=p)))
fits=c(fits,list(show.arima(x,order=c(p=5,d=3,q=4),p=p)))
